Day 7 Notes

Sampling and sampling bias

• Self-selection bias for surveys
• Problems with convenience samples: e.g., interview students at library on Friday night about their study habits
• Selection and survival bias: we don't see the kids who drop out of school, tracking a random sample of Alzheimer's patients to see how long they survive — you're more likely to sample those who survive a long time.

Random sampling

• Sampling frame and formal methods. Use sample

In-Class Activity

• Sampling Bias Activity s = select.books() which needs to be proceeded by

fetchData("simulate.r")

## Retrieving from
## http://www.mosaic-web.org/go/datasets/simulate.r

## [1] TRUE

Instructor's write up

Means and Models

Quick Review

• A model is a representation of something for a purpose.
• The statistical models that we're going to build in this class construct a function of zero or more inputs. The output of that function is the model value.
• By putting in the inputs for a given case, we get the fitted model value for that case.
• We fit models to a data set by building the function so that, for the data used for fitting, it gives a function that's as close to the data as possible. To give an idea about fitting …

fetchData("mLM.R")
mLM( wage ~ sex + sector, data=CPS85 )

• Each case has a residual, which is the difference between the actual value and the fitted model value. This represents “what remains unexplained” by our model.
• We can roughly judge the quality of a model by looking at the size of the residuals.
• The variance (function var()) provides a nice way to do this because it has a partitioning process.
  

Partitioning

The variance of the residuals plus the variance of the fitted equals the variance of the measured values.

This is much like the Pythagorean theorem.

It means we can measure “how much” variation has been “explained” by the model, and how much remains unexplained.

Try it on any grouping you like.

ACTIVITY: Construct a few models of the wage variable in the CPS85 data set, using different explanatory variables.

• Show that, in every one of your models, the variance partitions the wage into “explained” and “unexplained” parts.
• What did you conclude about wages from your various models?

Words for the Day

Accuracy versus Precision

• What's the difference?
• Synonyms for them.

The dictionary gets it substantially wrong:

precision |priˈsi zh ən| noun the quality, condition, or fact of being exact and accurate : the deal was planned and executed with military precision. • [as adj. ] marked by or adapted for accuracy and exactness : a precision instrument. • technical refinement in a measurement, calculation, or specification, esp. as represented by the number of digits given : this has brought an unprecedented degree of precision to the business of dating rocks | a precision of six decimal figures. Compare with accuracy . ORIGIN mid 18th cent.: from French précision or Latin praecisio(n-), from praecidere ‘cut off’ (see precise ).

Thesaurus synonyms: tools crafted with precision: exactness, exactitude, accuracy, correctness, preciseness; care, carefulness, meticulousness, scrupulousness, punctiliousness, methodicalness, rigor, rigorousness.

Today is all about Precision and how to measure and convey it. For us, the meaning of precision is “the degree of repeatability” or “reproducibility.” Sometimes this is conveyed by the “number of digits”, which is to say “the position of the last digit that matters”.

One wikipedia definition is useless: precision (statistics)http://en.wikipedia.org/wiki/Precision_(statistics).

Another one, accuracy and precision is more on target.

Review from last class

We saw that the choice of parameters in a model could be automated: fit the model, thereby choosing the parameters that eliminate bias in the residuals and minimize the variance.

Of course, the chosen parameters depend on the data to which the model is fitted. Insofar as this data is a random sample from a population, the parameters themselves are random. The question for use today is, “How random?”

The sampling process

Sample 100 runners from the runners data set and find the mean running time. Because there is some missing data, you'll need to do something like this:

options(na.rm=TRUE)

options(na.rm=TRUE)
run = fetchData("repeat-runners.csv")

## Retrieving from
## http://www.mosaic-web.org/go/datasets/repeat-runners.csv

nrow(run)

## [1] 24334

mysamp = sample(run,size=100)
mymod = mm( gun ~ sex, data=mysamp)
mymod

## 
## Groupwise Model Call:
## gun ~ sex
## 
## Coefficients:
##    F     M  
## 96.2  81.7

Across the class, comparing the different students, see what the range of values is.

Do it several times on the instructor's computer.

do(5) * mm( gun ~ sex, data=sample(run,size=100))

##        F     M sigma r.squared
## 1 101.31    NA    NA        NA
## 2  94.61 87.16 14.53    0.0629
## 3  97.95 86.06 12.70    0.1743
## 4  98.78 83.92 15.72    0.1840
## 5 100.27 84.55 14.02    0.2259

Do it many times with do and plot out the sampling distribution.

trials = do(100) * mm( gun ~ sex, data=sample(run,size=100))
densityplot( ~F, data=trials )

densityplot( ~M, data=trials )

Or, to put them both on one plot, which I want to do to compare to the population distribution, here's a more esoteric command:

require(reshape2)

## Loading required package: reshape2

densityplot( ~value, groups=variable, data=melt(trials[,1:2]))

## Using as id variables

How do we measure the width of the sampling distribution. With the standard deviation of the distribution.

sd(trials)

## Error: Error in is.data.frame(x): (list) object cannot be
## coerced to type 'double' Did you perhaps omit data= ?

This is called the standard error of the distribution.

IMPORTANT the standard error of the distribution is not at all the same as the distribution of the population. Here's that distribution of individual net running times:

densityplot( ~gun, groups=sex, data=run)

Resampling

You don't typically have the population at hand. (If you did … )

So we simulate the sampling distribution by constructing a simulated population that consists of many copies of the original. You could think of this as creating a pot of data, pouring into the pot many copies of the sample so that you have a huge set of cases to choose from.

In practice, this is implemented by sampling with replacement. To contrast sampling with resampling:

set = c(1,2,3,4)
sample(set)

## [1] 2 4 3 1

sample(set)

## [1] 4 2 1 3

sample(set)

## [1] 3 1 4 2

resample(set)

## [1] 1 2 1 1

resample(set)

## [1] 2 4 1 3

resample(set)

## [1] 1 2 4 4

The resample() function is set up, by default, to create a simulated sample that's the same size as the original.

Student task

Take your sample of size 100 and fit the model to a resampled set of data:

mm( gun ~ sex, data=resample(mysamp))

## 
## Groupwise Model Call:
## gun ~ sex
## 
## Coefficients:
##    F     M  
## 97.7  80.7

mm( gun ~ sex, data=resample(mysamp))

## 
## Groupwise Model Call:
## gun ~ sex
## 
## Coefficients:
##    F     M  
## 97.6  80.6

mm( gun ~ sex, data=resample(mysamp))

## 
## Groupwise Model Call:
## gun ~ sex
## 
## Coefficients:
##    F     M  
## 94.6  81.7

Do this many times to create the resampling distribution:

rs = do(100)* mm( gun ~ sex, data=resample(mysamp))

What's the width of this distribution?

sd(rs)

## Error: Error in is.data.frame(x): (list) object cannot be
## coerced to type 'double' Did you perhaps omit data= ?

That constitutes an estimate of the standard error of the coefficients of this model.

summary(mymod)

## Groupwise Model
## Call:  gun ~ sex 
## 
##   group 2.5 % 97.5 %
## 1     F 92.01 100.32
## 2     M 77.87  85.54
## 
## sigma:  14.2 
## r-squared: 0.208 
## Adj. r-squared: 0.2

QUESTION: Where does the information to do this come from, since mymod doesn't see anything but the particular sample that we took?

The Standard Error

The standard error depends on the size of the sample. More data means a more precise estimate. Let's quantify this.

Student Activity Use sampling and iteration for 200 times to calculate the standard error for different sample sizes. Assign a sample size to each student: 25, 50, 100, 200, 400, 800, 1600, 3200

Example:

trials = do(200)*mm( gun ~ sex, data=sample(run,size=800))
sd(trials)

## Error: Error in is.data.frame(x): (list) object cannot be
## coerced to type 'double' Did you perhaps omit data= ?

Draw a plot of the standard error versus sample size. Show that the standard error gets smaller as \( 1/\sqrt{n} \). You need 4 times as much data to get an estimate that's twice as precise.

From the Sampling/Resampling Distribution to a Confidence Interval

Two equivalent formats

• \( L \) to \( U \) with \( c \) confidence
• \( p \pm m \) with \( c \) confidence

Three components: 1. Point estimate \( p \) 2. Margin of error \( m \) 3. Confidence level \( c \) typically 95\%

Two main methods:

• Compute the 95\% coverage interval of the sampling distribution
• Make the margin of error 2 times the standard error. (Why 2? For bell shaped distributions it covers about 95\%.)

The confint program will do this whichever way you ask:

confint( trials, method="stderr")

##        name    lower   upper
## 1         F 95.73420 98.7618
## 2         M 84.59283 87.4524
## 3     sigma 14.19771 15.6386
## 4 r.squared  0.07703  0.1634

confint( trials, method="quantile")

##        name    2.5 %  97.5 %
## 1         F 95.59032 98.8736
## 2         M 84.68187 87.3467
## 3     sigma 14.30073 15.5605
## 4 r.squared  0.07492  0.1587

Simulation

Have students generate a random sample of 100,000 from a normal distribution with a specied mean and standard deviation. Then calculate the 95% coverage interval. Note that it extends from the mean about \( \pm 2 \) standard deviations.

Example:

samp = rnorm(100000, mean=20, sd=5)
qdata(c(0.025,0.975), samp) # should be about 10 to 30

##  2.5% 97.5% 
## 10.11 29.77

The second method, using the standard error, is very common, because it's easier to get a good estimate of the standard error than to get the 2.5 or 97.5 percentile.

Interpreting the Confidence Interval

Differences between groups

Simple visual method: do the confidence intervals overlap? Later on we'll see a more sophisticated method called hypothesis testing.

Coverage

Have everyone construct a groupwise model of the wage data against some grouping variable and pick some group to examine.

• First, do it to the entire data set, as if it were a population.
• Second, do it to a sample of 100 cases.
• Construct a 50% confidence interval.

Example:

popParameter = mm( gun ~ sex, data=run)
popParameter

## 
## Groupwise Model Call:
## gun ~ sex
## 
## Coefficients:
##  F   M  
## NA  NA

mymod = mm( gun ~ sex, data=sample(run,size=100))
confint(mymod,level=0.5)

##   group  25 %  75 %
## 1     F 96.37 99.27
## 2     M 87.15 89.82

Is the population parameter value within your confidence interval?

Have fun with the students whose interval doesn't cover their population parameter. Lead them to understand that with a confidence level of 50%, you'll cover the population parameter only half the time.

In any given study, you can't know if your confidence interval includes the population parameter. But if you follow the right techniques, your interval will include the population parameter 95% (or whatever the confidence level is) of the time.

Why not 100% confidence? This would have to extend to infinity.

Significant digits and the confidence interval

A simple method: Report two digits of the (doubled) standard error. Then report the point estimate to the same digit as the smallest digit in your standard error.

Examples:

• Point estimate=3418384.783, sd=67213.5.
  

Round doubled sd to two digits: 130000. Report point estimate to the same digit: 3420000 \( \pm \) 130000 with 95% confidence.

Going forward

Don't confusion precision with accuracy. The confidence interval tells you the precision of your estimate, based on the variability that comes about from the sampling process. But however precise your estimate is, it may still be inaccurate. Some reasons:

• Your sample may not be representative. Randomness fixes this, but not everything that you think is random is indeed random.
• There may be other factors involved that influence the system and that your model has not taken account of.
  

We'll be moving on to build models involving multiple explanatory variables in the next class. Part of doing that will be to expand the notion of models across groupwise means.

We saw an example of this in the first class, with the early studies of smoking and mortality. Breaking things up into groups by smoking status, without taking into account age, gave a result that we now know to be wrong.